﻿using System;
using System.Collections.Generic;
using Csu.Modsim.ModsimModel;
using ASquared;
using ASquared.ModelOptimization;
using ASquared.SymbolicMath;

namespace iwpm.optimization
{
    /// <summary>The base class that implements MODSIM's IConditionalGradientSolvableFunction.</summary>
    public class OptiFunctionBase : IConditionalGradientSolvableFunction
    {
        #region Instance variables

        private Model _model;
        private Link[] _links = null;
        private string[] _vars = null;
        private Symbol[] _symbolVars = null;
        private Symbol _fxn = null;
        private Matrix _mask = null;
        private double _scaleFactor;

        // decrease spills 
        private double _spillWeight = 1e4;
        private Link[] _spillLinks; 
        private Matrix _spillConstraint; 
        private LagrangeMultiplierMatrix _mu;
        private double _massBalSpillWeight = 5e8; 

        // decrease 
        private Link[] _routingLinks;
        private Link[] _routingLinks_Art;
        private Matrix _routingConstraint;
        private LagrangeMultiplierMatrix _lambda;

        #endregion
        #region Properties

        #region Linear problem feedback

        /// <summary>Gets and sets all the costs for all "decision links" (those links active in the optimization process) in the MODSIM _model.</summary>
        public Matrix costs
        {
            get
            {
                int num = _links.Length;
                double[] c = new double[num];
                for (int i = 0; i < num; i++)
                    c[i] = _links[i].mlInfo.cost;
                return c;
            }
            set
            {
                double[] c = value;
                int num = c.Length;
                for (int i = 0; i < num; i++)
                    if (c[i] > DefineConstants.COST_SUPERLARGE)
                        _links[i].mlInfo.cost = DefineConstants.COST_SUPERLARGE;
                    else if (c[i] < -DefineConstants.COST_SUPERLARGE)
                        _links[i].mlInfo.cost = -DefineConstants.COST_SUPERLARGE;
                    else
                        _links[i].mlInfo.cost = Convert.ToInt64(c[i]);
            }
        }
        /// <summary>Gets and sets the vector of decisions. In this case, decisions are flows in all the _links ordered according to Links_All.</summary>
        public Matrix decisions
        {
            get
            {
                int num = _links.Length;
                double[] q = new double[num];
                for (int i = 0; i < num; i++)
                    q[i] = _links[i].mlInfo.flow / _scaleFactor;
                return new Matrix(q);
            }
            set
            {
                if (value == null) return;
                for (int i = 0; i < _links.Length; i++)
                {
                    double flow = (value[i] * _scaleFactor);
                    if (flow > _links[i].mlInfo.hi)
                        _links[i].mlInfo.flow = _links[i].mlInfo.hi;
                    else if (flow < _links[i].mlInfo.lo)
                        _links[i].mlInfo.flow = _links[i].mlInfo.lo; 
                    else 
                        _links[i].mlInfo.flow = (long)flow;
                }
                _model.hydro.Update();
            }
        }
        /// <summary>Gets and sets the mask of nonlinear variables for the model network.</summary>
        public Matrix variableMask
        {
            get
            {
                return _mask;
            }
            set
            {
                _mask = value;
            }
        }
        /// <summary>Gets the vector of lower bounds. In this case, values are ordered according to Links_All.</summary>
        public Matrix lower
        {
            get
            {
                int num = _links.Length;
                double[] q = new double[num];
                for (int i = 0; i < num; i++)
                    q[i] = _links[i].mlInfo.lo / _scaleFactor;
                return new Matrix(q);
            }
        }
        /// <summary>Gets the vector of upper bounds. In this case, values are ordered according to Links_All.</summary>
        public Matrix upper
        {
            get
            {
                int num = _links.Length;
                double[] q = new double[num];
                for (int i = 0; i < num; i++)
                    q[i] = _links[i].mlInfo.hi / _scaleFactor;
                return new Matrix(q);
            }
        }
        /// <summary>Gets the variable names for this objective function.</summary>
        public string[] variables
        {
            get { return _vars; }
        }
        /// <summary>Gets the decision variables for this objective function in symbol form.</summary>
        public Symbol[] SymbolVariables
        {
            get { return _symbolVars; }
        }

        #endregion
        #region Nonlinear objective function

        /// <summary>A weight associated with the spill from the reservoir</summary>
        public double SpillWeight
        {
            get
            {
                return _spillWeight;
            }
            set
            {
                _spillWeight = value;
            }
        }
        /// <summary>A symbolic function of penalty terms to avoid spilling</summary>
        public Symbol SpillPenaltyFunction
        {
            get
            {
                Symbol s = 0.0, ms = 0.0;
                foreach (Link l in _links)
                    if (l.mlInfo.hydroSpillLinks != null)
                        for (int i = 0; i < l.mlInfo.hydroSpillLinks.Length; i++)
                        {
                            Link spill = l.mlInfo.hydroSpillLinks[i];
                            if (!spill.mlInfo.isArtificial)
                            {
                                s += (_symbolVars[spill.number - 1] ^ 2);
                                this.variableMask[spill.number - 1] = 1;
                            }
                        }
                foreach (Node n in _model.Nodes_Reservoirs)
                    if (n.mnInfo.spillLink != null)
                    {
                        int num = n.mnInfo.spillLink.number - 1; 
                        ms += _symbolVars[num] + (_symbolVars[num] ^ 2);
                        this.variableMask[num] = 1; 
                    }
                return _spillWeight * s + _massBalSpillWeight * ms;
            }
        }
        /// <summary>A parameter in the Lagrangian relaxation formulation of the routing alpha parameter</summary>
        public double RoutingAlphaParameter
        {
            get
            {
                return _lambda.Alpha; 
            }
            set
            {
                _lambda.Alpha = value; 
            }
        }
        /// <summary>A function that attempts to minimize the zero flow link to zero.</summary>
        public Symbol RoutingCompatibility
        {
            get
            {
                Symbol penalize = 0.0;
                _lambda.Update(_model.mInfo.Iteration, _routingConstraint);
                for (int i = 0; i < _routingLinks.Length; i++)
                {
                    Link r = _routingLinks[i], rArt = _routingLinks_Art[i];
                    Symbol q = _symbolVars[rArt.number - 1] - r.m.lagfactors[0] * r.mlInfo.flow;
                    if (_lambda.UseAugmented)
                        penalize += _lambda[i] * q + _lambda.Alpha / 2.0 * (q ^ 2);
                    else
                        penalize += _lambda[i] * q;
                    this.variableMask[rArt.number - 1] = 1;
                    _routingConstraint[i] = (double)rArt.mlInfo.flow 
                        - r.m.lagfactors[0] * (double)r.mlInfo.flow;
                }
                return penalize;
            }
        }
        /// <summary>The objective function of the MODSIM _model.</summary>
        public virtual Symbol fxn
        {
            get
            {
                if (_fxn == null)
                    throw new NotImplementedException("Need to set or override the objective function of the OptiFunctionBase class.");
                return _fxn;
            }
            set
            {
                _fxn = value;
            }
        }

        #endregion
        #region Convergence

        /// <summary>Allows setting of the MODSIM _model convergence fields during iterations of </summary>
        public bool IsConverged
        {
            get
            {
                return _model.mInfo.convgHydro;
            }
            set
            {
                _model.mInfo.convgHydro = value;
            }
        }

        #endregion

        #region Unnecessary

        /// <summary>This property is not needed. MODSIM serves as the linear optimization engine.</summary>
        public Matrix A
        {
            get { return null; }
        }
        /// <summary>This property is not needed. MODSIM serves as the linear optimization engine.</summary>
        public Matrix b
        {
            get { return null; }
        }
        /// <summary>This property is not needed. MODSIM serves as the linear optimization engine.</summary>
        public VariableType[] VariableTypes
        {
            get { return null; }
        }
        /// <summary>This property is not needed. MODSIM serves as the linear optimization engine.</summary>
        public ConstraintType[] ConstraintTypes
        {
            get { return null; }
        }

        #endregion

        #endregion

        #region Constructor

        /// <summary>The base objective function that is undefined.</summary>
        /// <param name="_model">The MODSIM _model.</param>
        /// <param name="parameters">The parameters for the particular objective function formulation.</param>
        public OptiFunctionBase(Model model)
        {
            _model = model;
            _scaleFactor = model.ScaleFactor;
        }

        #endregion

        #region Methods

        /// <summary>Initializes the instance for optimization.</summary>
        public virtual void Initialize()
        {
            _links = _model.Links_All; // Gets an array of links sorted by the link number
            _vars = new string[_links.Length];
            _symbolVars = new Symbol[_links.Length];
            for (int i = 0; i < _links.Length; i++)
            {
                _vars[i] = "q_" + _links[i].number;
                _symbolVars[i] = _vars[i];
            }
            _mask = new Matrix(_links.Length, 1.0);

            // build the routing compatibility function and links 
            List<Link> rLinks = new List<Link>();
            List<Link> rLinks_Art = new List<Link>();
            for (int i = 0; i < _links.Length; i++)
                if (_links[i].m.loss_coef >= 1) // check for routing links
                    for (LinkList ll = _links[i].to.OutflowLinks; ll != null; ll = ll.next)
                        if (!ll.link.mlInfo.isArtificial)
                        {
                            rLinks.Add(_links[i]);
                            rLinks_Art.Add(ll.link);
                        }
            _routingLinks = rLinks.ToArray();
            _routingLinks_Art = rLinks_Art.ToArray();
            _routingConstraint = new Matrix(_routingLinks_Art.Length, 1, 1.0);
            _lambda = new LagrangeMultiplierMatrix(_routingLinks_Art.Length);
            _lambda.UseAugmented = true;

        //    // build the spill link penalty
        //    List<Link> sLinks = new List<Link>();
        //    foreach (Link l in _links)
        //        if (l.mlInfo.hydroSpillLinks != null)
        //            for (int i = 0; i < l.mlInfo.hydroSpillLinks.Length; i++)
        //            {
        //                sLinks.Add(l.mlInfo.hydroSpillLinks[i]); 
        //                Link spill = ;
        //                s += (_symbolVars[i] ^ 2);
        //                this.variableMask[i] = 1;
        //            }
        //        return _spillWeight * s;
        //            private Link[] _spillLinks; 
        //private Matrix _spillConstraint; 
        //private LagrangeMultiplierMatrix _mu; 


        }
        public Symbol GetModsimObjectiveFunction(string[] excludeVariables)
        {
            // a list of variables to exclude in the objective function if desired.
            List<Symbol> exclude =
                excludeVariables == null ?
                new List<Symbol>() :
                new List<Symbol>(Array.ConvertAll(excludeVariables, s => (Symbol)s));

            // costs
            Matrix cost = this.costs;

            // get the function
            Symbol water = null;
            for (int i = 0; i < _symbolVars.Length; i++)
                if (!exclude.Contains(_symbolVars[i]))
                    if (water == null)
                        water = cost[i] * _symbolVars[i];
                    else
                        water += cost[i] * _symbolVars[i];
            return water;
        }

        #endregion
    }
}
